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Abstract 

In this paper, we apply Ritz variational approach to a new iso- 
chemical model of H2 molecule suggested by Santilli and Shillady. We 
studied Gaussian, Vg^ and exponential, 14, screened Coloumb poten- 
tial approximations, as well as the original Hulten potential, Vh, case. 
Both the Coloumb and exchange integrals have been calculated only 
for Ve. owing to Gegenbauer expansion while for Vg and Vh cases we 
achieved analytical results only for the Coloumb integrals. We con- 
clude that the 14-based model is capable to fit experimental data on 
H2 molecule in confirmation of the results of numerical HFR approach 
by Santilli and Shillady. Also, we achieved the energy-based estima- 
tion of the weight of the isoelectronium phase which is appeared to be 
of the order of 1%...6%, for the case of I4. However, we note that this 
is not the result corresponding to the original Santilli-Shillady model, 
which is based on the Hulten potential Vh- An interesting result is 
that in order to prevent divergency of the Coloumb integral for Vh the 
correlation length parameter Vc should run discrete set of values. We 
used this condition in our T4-based model. 







1 Introduction 



In this paper, we consider the four-body Santilh-Shillady isochemical model 
of H2 molecule [jl], |] characterized by additional short-range attractive 
Hulten potential between the electrons. This potential is assumed to lead 
to bound state of electrons called isoelectronium. The restricted three-body 
Santilli-Shillady model (stable and point-like isoelectronium) of H2 has been 
studied in ref . 0] , in terms of exact solution. For the mass of isoelectronium 
M = 2me, this solution implied much lower energy than the experimental 
one so we varied the mass and obtained that M = 0.308381me fits the experi- 
mental binding energy, Eexper[H2] = —1.174474 a.u. up to six decimal places, 
although at bigger value of the internuclear distance, R = 1.675828 a.u. in 
contrast to Rexper[H2] = 1.4011 a.u. We realize that the three-body model 
is capable to represent the binding energy but it is only some approxima- 
tion to the four-body model, and one should study the general four-body 
hamiltonian of the Santilli-Shillady model as well. 

In the present paper, we use Ritz variational approach to the four-body 
Santilli-Shillady isochemical model of H2 molecule, i.e. without restriction 
that the isoelectronium is stable and point-like particle, in order to find the 
ground state energy and bond length of the H2 molecule. 

In Sec. 2, we analyze some features of the four-body Santilli-Shillady 
isochemical model of H2 molecule. 

In Sec. 3, we apply Ritz variational approach to the four-body Santilli- 
Shillady model of H2 molecule. We calculate Coloumb integral for the 
cases of Hulten potential (Sec. 3.1.1), exponential screened Coloumb po- 
tential (Sec. 3.1.2), and Gaussian screened Coloumb potential (Sec. 3.1.3). 
Owing to Gegenbauer expansion, exchange integral has been calculated for 
the case of exponential screened potential, with some approximation made 
(Sec. 3.1.4). Exchange integrals for the Hulten potential and the Gaussian 
screened Coloumb potential have not been derived, and require more study. 
We present main details of calculations of the Coloumb and exchange inte- 
grals which have been appeared to be rather cumbersome, especially in the 
case of Hulten potential. 

In Sec. 3.2, we make numerical fitting of the variational energy for the 
case of exponential screened Coloumb potential Vg. Also, we estimate the 
weight of the isoelectronium phase. However, we use all the important results 
of the analysis made for the Hulten potential \4. 
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1) We conclude that the Ve-based model with the one-level isoelectronium 
is capable to fit the experimental data on H2 molecule (both the binding 
energy E and the bond length R). This is in confirmation of the results of 
numerical HFR approach (SASLOBE routine) to the V^-based model of ref. 

i- 

2) One of the interesting implications of the Ritz variational approach to 
the Hulten potential case is that the correlation length parameter r^, entering 
the Hulten potential, and, as a consequence, the variational energy, should 
run discrete set of values during the variation. In other words, this means that 
only some fixed values of the effective radius of the one-level isoelectronium 
are admitted, in the original Santilli-Shillady model, within the framework of 
the Ritz approach. This highly remarkable property is specific to the Hulten 
potential Vh while it is absent in the or V^-based models. 

3) Also, we achieved an estimation of the weight of the isoelectronium 
phase for the case of V^-based model which is appeared to be of the order 
of 1%...6%. This weight has been estimated from the energy contribution, 
related to the exponential screened potential Ve, in comparison to the con- 
tribution related to the Coloumb potential. 

4) Another general conclusion is that the effective radius of the isoelec- 
tronium Tc should be less that 0.25 a.u. 

We note that the weight of the phase does not mean directly a time share 
between the two regimes, i.e., 1...6% of time for the pure isoelectronium 
regime, and 99. ..94% of time for the decoupled electrons regime. This means 
instead relative contribution to the total energy provided by the potential 
and by the usual Coloumb potential between the electrons, respectively. As 
a consequence, the weight of the isoelectronium phase, which can be thought 
of as a measure of stability of the isoelectronium, may be 

1. Different from the obtained 1...6% when calculated for some other char- 
acteristics of the molecule, e.g., for a relative contribution of the pure 
isoelectronium to the total magnetic moment of the H2 molecule; 

2. Different from the obtained 1...6% for the case of the original Hulten 
potential Vh. 

So, the result of the calculation made in this paper is not the final result 
implied by the general four-body Santilli-Shillady model of H2 molecule since 
the latter model is based on the Hulten potential Vh- This paper can be 
viewed only as a preliminary study to it. However, we have made some 
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essential advance in analyzing the original Hulten potential case (Sec. 3.1.1), 
which we have used in the K-based model. 

Below, we describe the procedure used in Sec. 3 in a more detail. In 
Ritz variational approach, the main problem is to calculate analytically so 
called molecular integrals. The variational molecular energy, in which we are 
interested in, is expressed in terms of these integrals; see Eq.( p.2|) . These 
integrals arise when using some wave function basis (usually it is a simple 
hydrogen ground state wave functions) in the Schrodinger equation for the 
molecule. The main idea of the Ritz approach is to introduce parameters into 
the wave function, and vary them, together with the internuclear distance 
parameter i?, to achieve a minimum of the molecular energy. In the case 
under study, we have two parameters, 7 and p, where 7 enters hydrogen- 
like ground state wave function ( |3.10|) , and p = 'jR measures internuclear 
distance. These parameters should be varied (analytically or numerically) in 
the final analytical expression of the molecular energy, after the calculation 
is made for the associated molecular integrals. 

However, the four-body Santilli-Shillady model of H2 molecule suggests 
additional, Hulten potential interaction between the electrons. The Hulten 
potential contains two parameters, Vq and Vc, where Vq is a general factor, 
and is a correlation length parameter which can be viewed as an effective 
radius of the isoelectronium; see Eq. ( p.23| ). Thus, we have four parameters 
to be varied, 7, p, Vq, and Vc- The introducing of Hulten potential leads 
to modification of some molecular integrals, namely, of the Coloumb and 
exchange integrals; see Eqs. (|3.5|) and ( p.7|) . The other molecular integrals 
remain the same as in the case of usual model of H2, and we use the known 
analytical results for them. So, we should calculate the associated Coloumb 
and exchange integrals for the Hulten potential to get the variational energy 
analytically. In fact, calculating of these integrals, which are six-fold ones, 
constitutes the main problem here. Normally, Coloumb integral, which can 
be performed in bispherical coordinates, is much easier than the exchange 
one, which is performed in bishperoidal coordinates. 

Calculation of the Coloumb integral for Hulten potential, Vh, appeared 
to be rather nontrivial (Sec. 3.1.1). Namely, we used bispherical coordinates, 
and have faced several special functions, such as polylogarithmic function, 
Riemann (^-function, digamma function, and Lerch function, during the cal- 
culation. Despite the fact that we see no essential obstacles to calculate this 
six-fold integral, we stopped the calculation after fifth step because sixth (the 
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last) step assumes necessity to calculate it separately for each integer value of 
= {2'yrc)~^, together with the need to handle very big number of terms. 
During the calculations, we were forced to use the condition that should 
take integer values in order to prevent divergency of the Coloumb integral 
for Hulten potential. Namely, some combination of terms containing Lerch 
functions gives a finite value only if this condition holds. This condition is 
specific to Hulten potential. Note also that we can not get general form 
of a final expression for the Coloumb integral for Hulten potential because 
Lerch functions entering the intermediate expression (after the fifth step, see 
Eq.( |3.80| )) can be integrated over only for a concrete numerical value of their 
third argument. 

In order to proceed with the Santilli-Shillady approach, we invoke to two 
different simplified potentials, the exponential screened Coloumb potential, 
Ve, and the Gaussian screened Coloumb potential, Vg, instead of the Hulten 
potential Vh- They both mimic well Hulten potential at short and long 
range asymptotics, and each contains two parameters, for which we use the 
notation, A and Tc. In order to reproduce the short range asymptotics of 
Hulten potential the parameter A should have the value A = Vorc, for both 
the potentials. The Coloumb integrals for these two potentials have been 
calculated exactly (Sees. 3.1.2 and 3.1.3) owing to the fact that they are 
much simpler than the Hulten potential. Particularly, we note that the final 
expression of the Coloumb integral for Vg contains only one special function, 
the error function erf(z), while for Ve it contains no special functions at all. 

Having these results we turned next to the most hard part of work: the 
exchange integral. Usually, to calculate it one has to use bispheroidal coordi- 
nates, and needs in an expansion of the potential in some orthogonal polyno- 
mials, such as Legendre polynomials, in bispheroidal coordinates. Here, only 
the exponential screened potential Ve is known to have such an expansion but 
it is formulated, however, in terms of bispherical coordinates (the Gegenbauer 
expansion). Accordingly, we calculated exactly the exchange integral for Ve, 
at zero internuclear separation, i? = 0, at which case one can use bispherical 
coordinates. After that, we recovered partially the R dependence using the 
standard result for the exchange integral for Coloumb potential (Sugiura's 
result). Thus, we achieved some approximate expression of the exchange 
integral for the case of Ve- So, all the subsequent results correspond to the 
Ve-based model. 

Inserting obtained \4-based Coloumb and exchange integrals into the to- 
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tal molecular energy expression, we get the final analytical expression con- 
taining four parameters, 7, p. A, and Tc- Prior to going into details of the 
energy minimization for the T4-based (approximate) model, we analyze the 
set of parameters, and the conditions which we derived in the original Hulten 
potential case. 

(1) From the analysis of Hulten potential, wc sec (Sec. 2.1) that the 
existence of a bound state of two electrons (which is proper isoelectronium) 
leads to the following relationship between the parameters for the case of 
one energy level of the electron-electron system: Vq = h'^ /{2mr1). So, using 
the above mentioned relation A = VoVc we have A = 1/rc = 27/A, in atomic 
units {h = iTLe = c = 1). Thus note that, in this paper, we confined our 
consideration to the case of one-level isoelectronium. 

(2) Prom the analysis of the Coloumb integral for Hulten potential, we 
see (Sec. 3.1.1) that the condition, = integer number, should hold, and 
one can use it as well. 

We use the above two conditions, coming from the Hulten potential anal- 
ysis, in the energy minimization calculations for the case of our Vg-based 
model. The first condition diminishes the number of independent parame- 
ters by one (they become three, 7, p, and A) while the second condition means 
a discretization of the A parameter, A^^ = 4, 5, 6, . . . Here, we used the con- 
dition A~^ > 3 which we obtained during the calculation of the Coloumb 
integral for Vg. 

With the above set up, we minimized the total molecular energy of the 
Vg-based model. Numerical analysis shows that the A dependence does not 
reveal any minimum, in the interval of interest, 4 < A~^ < 60, while we have 
a minimum of the energy at some values of 7 and p. So, we calculated the 
energy minima for different values of A, in the interval of interest, 4 < A^^ < 
60. Results are presented in Tables 2 and 3. One can see that the binding 
energy decreases with the increase of the parameter r^ which corresponds to 
the effective radius of the isoelectronium. 

The following remarks are in order. 

(i) Note that the discrete character of Tc does not mean that the iso- 
electronium is some kind of a multilevel system, with different effective 
radia of isoelectronium assigned to the levels. We remind that we treat 

the isoelectronium as one-level system due to the above mentioned relation 
^0 = h^/{2mrl). In fact, this means that there is a set of one-level isoelec- 
tronia of different fixed effective radia from which we should select only one. 
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to fit the experimental data. 

(ii) The use of the exponential screened potential Ve can only be treated as 
some approximation to the original Hulten potential, and, thus, to the original 
Santilli-Shillady model of H2 molecule. So, the numerical results obtained 
in Sec. 3.2 are valid only within this approximation. Hulten potential makes 
a difference (one can see this, e.g., by comparing Sec. 3.1.1 and Sec. 3.1.2), 
and it is worth to be investigated more closely by, for example, combination 
of analytical and numerical methods. 

(iii) The results obtained in ref. [Q are based on the Gaussian screened 
Coloumb potential Vg approximation, to which the present work gives sup- 
port in the form of exact analytical calculation of the Coloumb integral for Vg 
(Sec. 3.1.3). Also, the present work gives possibility to make a comparative 
analysis of ref. [|l|, due to some similarity of the used potentials, Ve and Vg. 

(iv) Both the Coloumb integrals, for Vg and Vg, reveal a minimum in 
respect with A = 2'yrc, i.e. in respect with (see Figures 6 and 9) since min- 
imization in Ritz parameter 7 is made independently. In principle, this gives 
us an opportunity to minimize the total molecular energy Emoi with respect 
to Tc- However, there are two reasons that we can not provide this minimiza- 
tion. First, these minima correspond to rather large values of Vc, namely, 
Tc > 1 a.u. for Ve (Fig. 6), and Vc > 2 a.u. for Vg (Fig. 9). Of course, this 
is not an obstacle to do minimization but we note that we generally assume 
that the effective radius of the isoelectronium Tc is much less than the inter- 
nuclear distance, Vc R = Rexper[H2\ = 1.4011 a.u. Second, and the main, 
reason is that for the exponential screened potential case (Sec. 3.1.2) the pa- 
rameter A should be less than 1/3 to provide convergency of the associated 
Coloumb integral. Typically, 7 ~ 1.2, from which we obtain the condition 
Tc = A/27 < 0.2 a.u. Also, for the Hulten potential case (Sec. 3.1.1), we 
obtained A < 1/2, and hence Vc < 0.25 a.u. This means that, in fact, it is 
impossible to reach finite minimum of the total molecular energy E^oi in re- 
spect with Tc since the Coloumb integrals blow up, at > 0.25 a.u., leading 
thus to infinite total energy Emoi- So, in our approach we arrive at a strict 
theoretical conclusion that the effective radius of the isoelectronium rc must 
be less than 0.25 a.u. Clearly, this supports our assumption that Vc is much 
less than the internuclear distance R. 
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2 Santilli-Shillady model and the barrier 



In this Section, we consider the general four-body Santilh-Shillady model 
[H] of H2 molecule, in Born-Oppenheimer approximation (i.e. at fixed nu- 
clei). Shrodinger equation for H2 molecule with the additional short range 
attractive Hulten potential between the electrons is of the following form: 

-V? - TT-V^ - K,--— - + — (2.1) 



\ 2mi ^ 2m2 1 - e-'"i2A- 

ria rib R) 

where R is distance between the nuclei a and h. 

Interaction between the two electrons in the model is due to the potential 

V{n2) = Vc{ri2) + Vn{ri2) = — - K) ^ _ , (2.2) 

where ri2 is distance between the electrons, Vq and are real positive pa- 
rameters. Here, first term, Vc, is usual repulsive Coloumb potential, and the 
second term, V/^, is an attractive Hulten potential. 
Extrema of V{ri2) are defined by the equation 

V\ri2) = + = 0. (2.3) 

^ ' rl2 {e^-ih-- - 1)2 ^ ' 

In the limit ri2 00, potential V{rx2) ~ If\2 = Vc{ri2). Series expansion 
of V{ri2) at ri2 — > is 

V{n2)\r..^, = + I - ^n2 + O(rL). (2.4) 

In general, there is relationship of Hulten potential to Bernoulli polynomials 
Bn{x). Namely, Bernoulli polynomials are defined due to 

„^xs 00 



Y.Bn{x)-. (2.5) 



„=o 



and we can reproduce Hulten potential 



s 1 CO n 

-T^BnilA, (2.6) 
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taking s = —ri2/rc. First five Bernoulli coefficients are 

i?o(l) = 1, = i 52(1) = I, Bsil) = 0, B,il) = -1. (2.7) 

Eq. ( |2.6| ) means expansion of Hulten potential with the use of Bernoulli coef- 
ficients. 

Eq.( p.4|) implies that to have an attraction near vyi = 0, which is necessary 
for forming of isoelectronium, we should put the condition 

Vor, > e\ (2.8) 

We note that, in view of the asymptotics ( p.4| ), Q = ^/VqTc can be thought 
of as Hulten charge of the electrons. 

Under this condition, Viry}) has one maximum at the point defined by 
Eq. (|2.3|) . This is the equilibrium point at which the Coloumb potential is 
equal to the Hulten potential. So, we have barrier [B) separating two asymp- 
totic regions, (A) r — > and (C) r — oo, with Coloumb-like attraction and 
Coloumb-like repulsion, respectively. 

In the region A, attractive Hulten potential Vh dominates, and there- 
fore two electrons form bound state (isoelectronium), while in the region C 
Coloumb repulsion Vc dominates, and they are separated. This separation is 
limited by the size of the neutral molecule. For example, assuming that Hi 
molecule is in the ground state we have r < rmoi = 3.46 bohrs, where we have 
assumed that separation between the protons is R = 1.46 bohrs = 0.77ri2A. 

Existence of the bound state of the electrons and of the barrier B is a. 
novel feature provided by the model. The asymptotic states, in regions A 
and C, pertube each other due to the barrier effect in region B. 



2.1 Region A 

In the case 

Vor, > (2.9) 
we can ignore Coloumb repulsion Vc, and region A is a Hulten region, |V/i| 3> 



\Vc\; see Eq. (|2.4| ). Then, exact solution of one-particle Schrodinger equation 
with Hulten potential Vh, where wave function has the boundary conditions 
V^(0) = and ip{oo) = (see 0, problem 68), can be used to establish the 
relation between the parameters Vq and Tc, and to estimate Tc. 
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Energy spectrum for Hulten potential is given by 



?2 2 



= ' ^ = 1'2'---- (2-10) 



where 



f='-^rl (2.11) 



and m is mass of the particle. Assuming that there is only one energy level, 
namely, ground state n = 1, we obtain the condition 

= 1, (2.12) 

which can be rewritten as 

Note that this state is characterized by approximately zero energy, £^1 = 0, 
due to Eg. ( ^.101 ); strictly speaking, must be bigger but close to 1 in 
Eq. (|27[2|) . 

We should to note that the number of energy levels for Hulten potential 
is always finite due to Eq.( p.lO[) . Assumption that there are more than one 



energy levels in the bound state of two electrons, i.e. that /5 > 1, leads to 
drastical decrease of ground level energy Ei and relatively small increase 
of characteristic size of isoelectronium in the ground state. 

As the conclusion, the model implies "quantization" of the distance be- 
tween two electrons, r = ri2, namely, forming of relatively small quasiparticle 
(isoelectronium) characterized by total mass M = 2me, charge q = — 2e, spin 
zero, s = 0, and small size in the ground one-level state. This quasiparticle, 
as a strongly correlated system of two electrons, moves in the potential of 
two protons of H2 molecule, and one can apply methods developed for H2 
ion, with electron replaced by isoelectronium, to calculate approximate en- 
ergy spectrum of H2 0. However, this quasiparticle is not stable, being a 
quasi-stationary state, due to finite height and width of the barrier B. So, 
we must take into account effects of both regions B and C to obtain correct 
energy spectrum of H2 molecule, within the framework of the model. 
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2.2 Region B 



Quasiclassically, due to smooth shape of the barrier, and because of expo- 
nential decrease of wave functions inside the barrier, electrons are not much 
time in region so we can ignore this transient phase in subsequent consid- 
eration. 

We should to point out that the existence of the bound state in the region 
A and repulsion in the region C unavoidably leads to existence of the barrier. 

2.3 Region C 

In general, region C is infinite, rc < r < oo, where rc is the distance 
between two electrons at which the Hulten potential is much smaller than 
the Coloumb potential, \Vh\ ^ \Vc\- 

In this region, electrons are not strongly correlated, in comparison to 
that in region A. Here, correlation is due to usual overlapping, Coloumb 
repulsion, exchange effects, and Coloumb attraction to protons. Shortly, we 
have the usual set up as it for the standard model of H2 molecule. 

Discarding, for a moment, effects coming from the consideration of regions 
A and B, we have finite motion of the electrons in region C. Namely, in the 
ground state of H21 the distance between electrons is confined by r = Vmoi — 
3.46 bohrs. We restrict consideration by the ground state of H2 molecule. 

Due to this finiteness of the region C, r < r„ioU two electrons on the same 
orbit have constant probability to penetrate the barrier to form strongly 
correlated system, isoelectronium, and vice verca. 

2.4 Model of decay of isoelectronium 

Below, we assume that the isoelectronium undergoes decay, and the resulting 
two electrons are separated by sufficiently large distance, in the final state. 
This leads us to consideration of the effective life-time of isoelectronium. To 
estimate the order of the life-time, we use ordinary formula for radioactive a- 
decay since the potential V{r) is of the same shape, with very sharp decrease 
at r < r^nax and Coloumb repulsion at r > r^ax- This quisiclassical model is a 
crude approximation because in fact the electrons do not leave the molecule. 
Moreover, we have the two asymptotic regimes simultaneously, with some 
distribution of probability, and it would be more justified here to say on 
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frequency of the decay-formation process. However, due to our assumption 
of small size of isoelectronium, in comparison to the molecule size, we can 
study an elementary process of decay separately, and use the notion of life- 
time. 

Decay constant is 




(2.14) 



where we put, in atomic units, 

h^l, e = 1, m = 1/2, Z = 1, rmax = 0.048, E ^ 1. (2.15) 

Here, E = 1 a.u. = 27.212 eV is double kinetic energy of the electron on first 
Bohr's orbit, Qq — 0.529ri2^, that corresponds approximately to maximal 
relative kinetic energy of two electrons in ground state of H2, and m — 1/2 
is reduced mass of two electrons. 

We obtain the following numerical estimation for the life-time of isoelec- 
tronium: 

l/A = Do-1.6-10-^^sec, (2.16) 

i.e. it is of the order of 1 atomic unit of time, r = 2.42 • 10"^'' sec. For lower 
values of the relative energy E, we obtain longer lifetimes; see Table 2. 

The quasiclassical model for decay we are using here is the following. 
Particle of reduced mass m = 1/2 penetrate the barrier B. This means a 
decay of isoelectronium. In the center of mass of electrons system, electrons 
undergo Coloumb repulsion and move in opposite directions receiving equal 
speed so that at large distances, r ^ Tmaxi each of them have some kinetic 
energy. This energy can be given approximate upper estimation using linear 
velocity of electron on first Bohr's orbit, v = 2.19-10^ cm/sec, since electrons 
are in the ground level of H2 molecule (this is the effect of the nuclei). This 
upper estimation corresponds to assumption of zero velocity of the center 
of mass in respect to protons which we adopt here. Kinetic energy of the 
particle of reduced mass is then double kinetic energy of electron, in center 
of mass system. 

As the conclusion, in the framework of the model, H2 molecule can be 
viewed as a mixed state of H2 ion like system, i.e. strongly correlated phase 
(Hulten phase), when electrons form isoelectronium, and standard model of 
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Energy E, a.u. 


eV 


Lifetime, Do-sec 


2 


54.4 


2.6-10-1^ 


1 


27.2 


1.6 ■ 10-1^ 


0.5 


13.6 


2.2 ■ 10-1** 


0.037 


1 


5.1 • 10-6 


0.018 


0.5 


4.0 


0.0018 


0.1 


3.1 • 10+25 



Table 1: Lifetime of isoelectronium. E is relative kinetic energy of the elec- 
trons, at large distances, r ^ Tmax, in the center of mass system. 



H2, i.e. weakly correlated phase (Coloumb phase), when electrons arc sepa- 
rated by large distance, r > rmax- Note that, as it has been mentioned above, 
we ignore the transient phase (inside the barrier) in this consideration. Ev- 
idently, the (statistical) weight of each phase depends on the characteristics 
of the potential V{ri2). 

For extremally high barrier, only one of the phases could be realized with 
some energy spectra in each phase, namely, either spectrum of H2 ion like 
system (with electron replaced by isoelectronium), or usual spectrum of H2 
molecule (without Hulten potential), respectively. 

For high but finite barrier, each phase receives perturbation, and their 
(ground) energy levels split to two levels corresponding to simultaneous re- 
alization of both the phases. Note that the value Vmax is indeed high, 
Vmax ^ 500 eV, under the given values of the parameters. 

In general, existence of the strongly correlated phase (isoelectronium) 
leads to increase of the predicted dissociation energy, D, of H2 molecule. 
Indeed, the mutual infucncc of the regions A and C decreases the ground 
energy level E of H2 due to the above mentioned splitting. The general 
formula for D is 

D = 2Eo- {E + -huj), (2.17) 

where 2£'o = — 1 is total energy of two separated H atoms, and )^hio is zero 
mode energy of the protons in ■ So, decreasing of £J < causes increase 
oiD. 

It is remarkable to note that experimental data give dissociation energy 
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Dexper[H2\ = 4.45 eV for H2 molecule (see, e.g. and references therein) 
while theoretical predictions within the standard model are D = 2.90 eV 
(Heitler-London), D = 3.75 eV (Flugge), and D = 4.37 eV (Hylleraas). We 
observe that improvement of the variational approximation gives better fits 
but still it gives lower values (about 2% lower) partially due to the fact that 
variational technique used there predicts generally bigger value (upper limit) 
for the ground energy. 

Below, we use the same Ritz variational technique as it had been used by 
Heitler, London and Hylleraas but the feature of the model is the existence of 
additional attractive short range potential between the electrons suggested 
by Santilli and Shillady. 

3 Variational solution for ground state en- 
ergy of H2 molecule 

In the limiting case of large distances between the nuclei, R — > 00, we have 
the total wave function of the electrons in the form 



where the first term corresponds to the case when electron 1 is placed close to 
nucleus a and /(r^i) is wave function of the corresponding separate H atom 
while the second term corresponds to the case when electron 1 is placed close 
to nucleus b. Symmetrized combination ('+' sign) corresponds to antiparal- 
lel spins of the electrons 1 and 2, and, as the result of the usual analysis, 
leads to attraction between the H atoms. Below, we use this symmetrized 
representation of the total wave function as the approximation to exact wave 
function. 

3.1 Analytical calculations 

By using Ritz variational approach and representation ( p.l| ) , we obtain from 
the Schrodinger equation ( |2.1| ) the energy of H2 molecule in the following 




(3.1) 



form (cf. §), 



A + A'S _ 2{C + £S) - jC' + S') 1_ 
1 + 52 1 + 52 ^ R' 



(3.2) 
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where 



is overlap integral, 



S = Jdv nrai)f{ni) (3.3) 



C= ldv—\f{ra,)\\ (3.4) 



// p-'Vlj'c \ 

dvidv, [—-VoY^^z^j\f{rai)\Vin2)\', (3.5) 

are Coloumb integrals, 



S= I dv —nrai)f{n,), (3.6) 



E' = j dv^dv^ f — - V^^—^^\ /*(r,i)/(rw)/*(r.2)/(r,2) (3.7) 
are exchange integrals, 

A = \dv rir,,) (-^V? - i-) /(r,i) (3.8) 

and 

A' = ldv /*(r,i) (-iv? - ^) /(r,0. (3.9) 

We use atomic units, e = 1, mi = m2 = me = 1. 

Quite natural choice is that the wave functions in Eg. ( p.l|) are taken in 
the form of hydrogen ground state wave function, 

f^r) = ^e-^^ (3.10) 

where 7 is Ritz variational parameter (7=1 for the proper hydrogen wave 
function), and r = r^i, r{,i, ra2, Tfe2- With the help of 7 we should make better 
approximation to an exact wave function of the ground state. Namely, we 
should calculate all the integrals presented above analytically, and then vary 
the parameters 7 away from the value 7 = 1 and R in some appropriate 
region, say 1 < i? < 2, to minimize the energy (3^). As the energy min- 
imum will be identified the found value of the parameter R corresponds to 
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1 2 3 4 5 6 

r , a . u . 



Figure 1: The overlap integral 5 as a function of p, Eq. Here, p = ■jR, 

where 7 is Ritz parameter and R is the internuclear distance. 

optimal distance between the nuclei. This value should be compared to the 
experimental value of R. 

All the molecular integrals (p.3|) - (|3.9|) , except for the Hulten potential 
parts in ( P75| ) and ( |3.7| ), are wellknown and can be calculated exactly; see, 
e.g. 1^. Namely, they are 

^5 = (1 + P + ^p') e-", (3.11) 

C^Cc = ^il-{l + p)e-'n, (3.12) 
P 

C'c - C|Vo=o = - (1 + + i/^' + Ip'^^'"")^ (3.13) 

S = 80 = 1(1 + p)e-', (3.14) 

h{p) = S\p){\np + C) - S'{-p)E,{4p) + 2S{p)S{~p)E,{2p), (3.16) 
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Figure 4: The exchange integral S'f^ as a function of p, Eq. ( p.l5| ). 



e ^ 



E^ip) = I —dt, (3.17) 



p 



^=-^7' + 7(7-1), ^'=-lf5 + ^(^_i^^^ ^3_18) 
where C is Euler constant, and we have denoted 

p = 7i?, (3.19) 

which can be taken as a second Ritz variational parameter in addition to 
7. The most hard part of work here is the exchange integral (|3.15|) , which 
was calculated for the first time by Sugiura (1927), and contains one special 
function, the exponential integral function Ei{p). 

Our problem is thus to calculate analytically the Hulten potential parts 
of the Coloumb integral ( p.5| ) and of the exchange integral ( p.7| ), and then 
vary all the Ritz variational parameters in order to minimize the ground state 
energy (jS^), 

-Emoi (P^^^^sters) = minimum. (3.20) 

In general, we have four parameters in our problem, Emoi = Emoiil, p, Vq, Vc), 
with the first two parameters characterizing inverse radius of electronic orbit 
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and the internuclear distance, respectively, and the last two parameters com- 
ing from the Hulten potential. However, assuming that the isoelectronium 
is characterized by one energy level, i.e. /3 = 1, we have the relation (^]T^) 



between Vq and so that we are left with three independent parameters, say, 
p,rc). In fact, we have three independent parameters for any 
fixed number jS of the levels due to the general relation ( p.llj ). 

Behavior of the energy E^oi as a function of 7 and p is more or less clear 
owing to known variational analysis of the standard model of H2 molecule. 
Namely, Emoi reveals a local minimum at some values of 7 and p. Thus, we 
should closely analyze the Vc dependence of the energy which is specific to 
the Santilli-Shillady model of H2 molecule. 

Below, we turn to the Coloumb integral for the Hulten potential. 

3.1.1 Coloumb integral for Hulten potential 

To calculate the Hulten part of the Coloumb integral ( p.5|) we use spher- 
ical coordinates, {rf,2,92,(p2), when integrating over second electron, and 
{rbi,9i,{pi), when integrating over first electron. 
The integral is 



7r 00 7r 00 



C'j, = 4n^J d9, Jdrbi Jd92 Jdn2 Vh{n2){^e-^^''''n2' sm92)x (3.22) 



0000 



V 71- / ' 

where Hulten potential is 



Vh{n2) = , , , =-■ (3.23) 



Here, we have used 



rai = ^Jni^ + R"^ -Hi'^R cos 9i , 
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ri2 yr-fea^ + rfci^ - 2rb2rbi cos6'2, 

and the fact that integrals over azimuthal angles Lpi and give us 47r^. 
First, we integrate over coordinates of second electron, 

TT OO 2 

1 = 271 J J dn2 Vhiru) (^e-^^'-^^na' sin 62) . (3.24) 


Integration over 62 gives us 

1^1 dn2 (h + l2 + l3 + h + h), (3.25) 



where 



and 



7l = _4^3g-27r,2^^22^ (326) 

72 = -27V,— ./(rbi - r;,2)2e-2^''''2 ln(l - eV('-w-'^02)'AcN (3 27) 

73 = -27Ve— A/(r6i + n2ye-^^"'' ln(l - eV(^^^+^/'^=), (3.28) 

74 = 27 "-^e-^^'"' Li2 (eVi^bi-n2)Vr,\ (339) 

h = 27Ve'— e-'^''''2Li2(eV('^w+'-''2)V'^cN (330) 

00 ^fe ° 



^^2(.) = E S = / i^^^^^ (3.31) 

is dilogarithm function. 

Now, we turn to integrating over rb2. For 7i we have 

00 

jdn2h^-l. (3.32) 


In 72, we should keep (r^i —rb2) to be positive so we write down two separate 
terms, 

00 Tfji CO 

j drb2h = hi + l22 = j drb2 I2{rb2 < Hi) + j drb2 h{rb2 > ni). (3.33) 
n-i 



19 



In these two integrals, I21 and I22, we change variable to x and y, respec- 
tively, 

X = {m - rb2)/rc, ni/vc < x < 0, y = {n2 - ni)/rc, < y < cx), 

(3.34) 

in order to simplify integrating. In terms of these variables, we have 


721 = I dx 2^^rMx - JjLx^)e~^iirbi-rcx) ^^f^^ _ gX)^ ^3 35) 

J Th-l 



Hi Ac 



00 

722 = - / dy 2'y\c\y + ^|/2^)g-27(ra+rcj/) ^^^^^ _ ^yy (3_3g) 







We are unable to perform these integrals directly. To calculate these inte- 
grals we use method of differentiating in parameter. Namely, we use simpler 
integrals, 

Li^ Jdx e'^'^^^ ln(l - e^) (3.37) 

and 

L2 = jdy ln(l - e^), (3.38) 

and differentiate them in parameter Tc to reproduce I21 and 722 ■ (One can 
use parameter 7 for this purpose, or introduce an independent parameter 
putting it to one after making calculations, with the same result.) Namely, 
by using definitions of Li and L2 we have 

- '^''■Ki:'' - 4^|^^'> 1=°../- p-^^) 



'2 



Now, the problem is to calculate indefinite integrals, Li and L2, which make 
basis for further algebraic calculations. After making the calculations, we 
have 

^1 = 2-^e'"''^(27rc($(e^ 1, 2^r,) + ln(l - e^)) - l) (3.41) 
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and 

^2 = -^^e-2^^^(27re(<^>(e^ 1, -2^r,) + ln(l - e^)) + l), (3.42) 



where 



^{z,s,a) = Y. a + fc^O, (3.43) 

fc=o + 

is Lerch function, which is a generahzation of polylogarithm function Lin{z) 
and Riemann (^-function. Particularly, Li2{z) = $(2;, 2,0). Also, we note 
that the Lerch function arises when dealing with Fermi-Dirac distribution, 
e.g., 

°° us 

I gfe-M + i = + s + 1, 1). (3.44) 



Below, we will need in derivatives of Lerch function $(2;, s, a) in third argu- 
ment. By using the definition ( |3.43|) we obtain directly 

-^$(2;, s, a) = $'(z, s, a) = -s^(z, s + 1, a), (3.45) 
da 

(P 

— $(^, s, a) = $"(2, s, a) = s{s + 1)<I'(^, s + 2, a). (3.46) 
Inserting (|^) and (g^ into ( 1339| ) and (|^ we get 



^ (6-27(^.1 -r-ex) (3 ^ 27(r,i - (2 + 7rw)r,x + Tr^'a^^)- (3.47) 

■vrLi V 

-27r,((l + 7(rfei - 2(1 + ^ni)r,x + 27r,2x2)) [$(e", 1, 27re) + ln(l - e")] + 



hi- . 

47^^61 



+27r,(-(l + 7(r,i - 2r,x))<l''(e^ 1, 27rJ + ^r,V{e^, 1, 27rJ)) 



ix-O 



/22 = -^(e-2^(^"+^^^)(3 + 27(rw + (2 + ^n^)r,x + ^r,\^)+ (3.48) 

+27re((l + 7(r6i + 2(1 + -ini)r,x + 27r,2x2))[<^)(e^ 1, -27r,) + ln(l - e")] + 

+27r,((l + 7(r,i + 2r,x))$'(e^ 1, -27r,) + ir,^'\e\ 1, -27rJ))) \lz^ . 

Now, we have to use the above derivatives (|3.45|) and (|3.46|) of Lerch function 
to obtain final expressions for I21 and l22- Then, we should take the limits 
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X — >■ rfei/rc, X — > 0, and ?/ — 0, ?/ — oo, respectively. The endpoints 
X = Tfji/rc and y = oo can be inserted easily, with the endpoint y = oo 
yielding zero, while the limits x — and y ^ require some care because 
of the presence of some divergent terms. 

To collect all the terms, we sum up J21 and (— 1)/22 given by ( p. 47] ) and 
(|3.48|) , put X = y, and take common limit a; — > 0, inserting x = for 
polynomial and exponential (welldefined) terms. We get 

= --^(ree-2^'-"(27re(l + 7rw)[$(e^2,27r,) - $(e^ 2, -27re)]+ (3.49) 

+47V[$(e", 3, 27r,) + $(e", 3, -27r,)] + 

+ (1 + 7^61) [$(e^ 1, 27r,) + <l>(e^ 1, -27r,) - 2 ln(l - e")]) 

The limits of Lerch functions of second, $(e^, 2, ±27rc), and third, 

$(e^', 3, ±2'jrc), order, at x — >■ 0, are welldefined while each of the terms in 

S(27rJ = lim[$(e'^, 1, 27rJ + ^{e"", 1, -27rc) - 2 ln(l - e"^)] (3.50) 

is divergent since Lerch function of first order, $(e^', 1, ±27rc), increases un- 
boundedly at X — >■ 0. We will analyze this limit below, to identify the con- 
dition at which the divergencies cancel each other. Now, we collect all the 
terms obtaining final result for the integral in the form 

[dri,2 h = ^(-(-3 + 27rbl + 87V$(e'■"/'•^3,27r,)+ (3.51) 

i Am ^7 

+27rc(l - 7r,)[$(e'^"/'^^ 1, 27r,) + 2-fr^(l>{e'''^^^''% 2, 27re) + ln(l - e''"/''^)])- 
-2ree-2^^" (1 + 7rw){5(27r,) + 27r,[C(2, 27re) - C(2, -27re)] + 
+47V[C(3, 27re) + C(3, -2ir,)]}), 

where 

k=l + f^) 
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is generalized Riemann (^-function. The values of ({2, ±2'yrc) and ({3, ±2'yrc) 
entering ( p.51|) are welldefined. For example, at 7 = 1.4 and Tc = 0.0048, we 
have 

C(2, ±27rJ ~ 5537, C(3, ±27re) - 2462. (3.53) 

Now, we turn to close consideration of the limit ( |3.50| ) entering (|3.51| ). 
Let us calculate it for the particular value 27rc = 1/100. Using expansion of 
each term of B around a; = 0, we obtain 



1 



101. _ 1 



B{—) = lim 100 - ^^{100Ti—)iC + ln(l - + ^(— ))}- (3.54) 



r(^)^ ^100^^ ^ ' ^^100' 



1 199 99 

-^{100r(— )(C + ln(l -s)+ ^(— ))} + 2 ln(l -s)+0{l~s] 



99r(^) ' '100'' ' ' "^'100 



X 



where we have denoted, for brevity, s = e 

°° 1 r'(z) 

m = E ^ = WT (3.55) 

is digamma function, T{z) is Euler gamma function, and C is Euler constant. 
Using elementary properties of gamma function we obtain from Eq.( |3.54]) 

1 1 99 
B( ) = 100-2C-^( )-^(— ), (3.56) 

so one can see that the logarithmic divergent terms cancel each other, and 
the limit is welldefined for 27rc = 1/100. The same is true for any integer 
value of 

k = J- (3.57) 

27rc 

while at noninteger k the limit -B(^) blows up. Generalizing the above par- 
ticular result ( p.56| ), we can write down 

5(1) = k-2C- ^(^) - ^(1 - ^), (3.58) 

for any integer k > 2. 

This highly remarkable result means that to have finite value of the 
Coloumb integral we should use the condition that A^^ = (27rc)^^ = k 
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is an integer number. Recalling that typically 7 ~ 1.5 and Tc — 0.01 we have 
the integer number k ~ 30. 

Now, we turn to the next integral, I3. It is similar to I2 so that we present 
the final expression, 

CO 

/ drb2 h = 1 (3+27r;.i+27r,(l+7rw)[$(e'^"/'^=, 1, -27re)+ln(l-e'-"/'^=)] 

(3.59) 

-47Vc^(l + 7rw)$(e^"/''% 2, -27rc) + 87Vc^$(e''"/''% 3, -27rc)). 
The integral 1^ is more complicated, 

00 

j drb2 h = hi + ^42, (3.60) 



where 



I dn22jW—e-^^'''Li2{e^''^'-'^^^'> (3.61) 
i '''' 
00 

742= f dn22-fW—e-^^''''Li2(e^'-''-''''^H. (3.62) 

Introducing variables 

X = (rw - rb2)/rc, y = {^2 - ni)/rc, (3.63) 

we rewrite the integrals in the form 


741= / c^x 27Ve'^^'"""'"=''^[Li2(en - — xLi2(e^)], (3.64) 

Tbi/rc 
00 

/42 = - / dy 27Ve'e2^('^«+'-=^)[Li2(e^) + — yLi2(e^)]. (3.65) 

i 

In the r.h.s. of /41, the first term can be calculated directly in terms of 
Lerch function while the second term can be obtained from the first term by 
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differentiating it in the parameter, for which we choose again Tc- Namely, 
the basic integral, which we will use to calculate /41, is 



Ml = Jdx e2^"="L^2(e"), (3.66) 

XO 

for which we have 

+r(l + 2-fr,){jry - 3C - 3V'(27r,))- 
-3e2™r(l + 27r,)(<l>(e^«,l,27r,) +ln(l -e^«) + 27r,Li2(e^''))). 



We use this result in the first term of /41. Differentiating Mi given by Eqs. 
(|3.66|) and ( p.67|) in Vc, we reproduce the second term of /41, up to a factor. 
So, collecting these results and inserting xq = r^i/rc we obtain after some 
algebra 

hi = - ^2n,Til + 2^r,) (""'""["-^9 + 4vrVrcV,i)r(27re)+ (3.68) 

r(l + 27rc)(6Crc - 3rbi - 6C7rcr;,i - vr^r^^ - 6(7rbi - l)rc^(27rc)- 

-67r,V(27'^c))] + 3(2rwr(l + 27r,) - 3rer(27r,) + 

+2rer(l + 27r,)[(l - 2^ni)^e'>''/''% 1, 27r,) + 7r,$(e'-"/^% 2, 27re) + 

+ (1 - 27^1) ln(l - e^'"/'-^) +7(1 -47r,i)r,7/;'(e^'"/^'^)]), 

where ip'{z) = dip{z)/dz is derivative of digamma function. 

To calculate I42 we use a similar method. However, care should be exerted 
when taking limit ?/ — 0. The basic integral, which we will use to calculate 
^42, is 

M2 = - J dy e^""^yLi2{e-y), (3.69) 

where we have replaced y —y so that the endpoints will be due to < 
y < —00. The result for M2 is 

M2 = -^^e'^^"'^y(l + 27r,e^$(e^ 1, 1 + 27r,) + 27r,ln(l - e~^)- (3.70) 
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We should insert here the endpoints ?/ = and y = — oo. In the hmit 
y — > —oo, M2 is zero. In the limit y —>■ 0, we have 

Lt,ie-y)\y^o = Y (3-71) 
and, assuming that k = 1/ {2'yrc) is an integer number, 

$(e^ 1, 1 + 27r,) + ln(l - e-^)|,_o = -(-^ + C + tlj{2^r,)). (3.72) 



Thus, 



Mal^Io °° = Y^^i^C + vrVc + 3^(27r,)), (3.73) 



for integer k. We should point out that, in the case of noninteger k, M2 
increases unboundedly at y ^ 0. Using this result in J42, we obtain 

/42 = -Ti^e-27-« (6C(1 + 7rw) + vr Vc + 2n^-f\,ni+ (3.74) 



+6(1 + 7rbi)?/'(27rc) - 67rc^'(27rc)^ 

Summing up /41 given by ( p.68|) and I42 given by (|3.74|) , we get 

00 

/ cirb2 h == --— — — — (3r,r(27rJ + e-2^^i-r,ir(l + 27r,)+ (3.75) 







+r,(-3r(27re) + 4r(l + 27re)(-(l + ^ni)iC + ^(27r,) + 7r,^'(27r,)))]- 
-2rj:{l + 27rJ[<l>(e''"/''% 1, 27rJ + ln(l - e'^"/''^ + 7rc($(e''"/'^% 2, 27rJ + 

+L^2(e'-"/'-=))]). 

The integral Is is similar to I4 so that we present the final expression, 



[dn2 h = -^^(3 + 47re[e-^"/"^($(e-""/"% 1, 1 + 27re)+ (3.76) 
+$(e-''"/''% 2, 1 + 27r,)) + ln(l - e""/'^) - ^r,ip'{e""/''^)]) . 
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Now, we are in a position to sum up all the calculated integrals Ji, . . . , J5, 
and obtain, due to ( p.25 ), the following final expression for the Coloumb 
integral over coordinates of second electron. 



Hm) = -( 



7l(l + 



)ctg(27rc7r)e' 



2 87rw 
1 



7^61 



mi 



)e-2^""+ (3.77) 

e-'^i-i/rc^^g-rw/r,^ 1, 1 + 27re) + 

+ 



I 2 2 

+7 Tc 



e-r-w/rc^^g-rw/r-c^ 2, -27r,) - <l>(e""/"% 2, -27r,)- 



1 ,1 



-<l>(e""/"% 2, 27re) + (-<l>(e""/"% 2, 27r,) - $(e''"/'^% 2, -27re)) + 

7^61 2 

+^e-27'-.i^'(27r,) + 6^27^.1 + ^)(C(2, -27r,) - C(2, 27r,)) 
mi 7rfei 



+ 



+ 



7^61 



<l>(e""/"% 3, 27re) + $(e""/"% 3, -27r,)- 



-^^^•"(C(3,27r,) + C(3,-27r,)) 



where we have collected the terms due to power degrees of r^. It should be 
stressed that here (27rc)~^ is assumed to be an integer number. The above 
expression represents the Hulten part of the electrostatic potential caused by 
charge distribution of the second electron. 

Next step is to integrate ( |3.77| ) over the coordinates of first electron. 



7r 00 




Prior to that, we denote 



A = 27rc = -, r = 7^1, 



(3.79) 
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and rewrite Eq. ( p.77| ) in a more compact form, 



/(r) = -il + |-)e-2'- + \x\7r{l + -)ctg(7rA)e-2^' - le-2'-A$(e~2^'/\ 1, 1 + A) + 

(3.80) 



+$(e2^/\ 1, -A) - <l>{e^'/\ 1, A) + -$(e2^/\ 1, -A) 

r 



+ 



+— e-2^'/^$(e-2'^/\ 2, -A) - <l>(e'^/\ 2, -A) - $(e''"/\ 2, A) + 
4 L 2r 

-(^$(e2'^/\ 2, A)-<l>(e2^/\ 2, -A)+e-2>'(A))+e-2^(l+-)(C(2, -A)-C(2, A))' 



+ 



4r 



$(e2^/\ 3, A) + ^{e^^'\ 3, -A) - e-2^(C(3, A) + C(3, -A)) 



Since /(r) does not depend on 9i one can easily integrate over 9i in Eq. 
and then change variable r^i to r = 7r;,i, obtaining 



C'^ = — I dr I{r) [(1 + 2Y/(p-r)2)re-V(p-'-)' - (1 + 2(p + r))re-2(''+") 

(3.81) 

where p = •jR. Again, we should use separate intervals to keep (p — r) to be 
positive, namely, we rewrite as 



where 



C'l^ — Ji + J2 + J3, 

1 ? 

= — dr /(r)(l + 2p - 2r)re" 



-2(p-r) 



J2 



— y"c/r /(r)(l + 2r-2p)re-^("-''). 



1 °° 

J3 = -— J dr /(r)(l + 2p + 2r)re~ 



2{p+r) 



(3.82) 

(3.83) 
(3.84) 
(3.85) 



Now, we are ready to make integration over the last remaining variable, r, to 
obtain complete analytical expression of the Coloumb integral for the Hulten 
potential. 
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However, Lerch functions entering Eq.( p.80D make obstacle to do integral 
(|3.82 ) for a general case because they have different functional form for dif- 



ferent values of the parameter A. So, each of the above integrals Ji,2,3 should 
be calculated independently for every numerical value of A. Moreover, for 
the values of interest, e.g., A = 1/30, each Lerch function is expressed in the 
form of sum of elementary functions with too big number of nontrivial terms 
to handle them (incomplete Euler beta function arises here). So, the integral 
cannot be reliably calculated even for a single value of A, within the interval 
of interest, A = 1/30, 1/31, . . . , 1/100. Also, elementary analysis shows that 
we can not implement the assumption of small Vc into Eq.( p.80| ), to use first 
order approximation in Tc. Indeed, Lerch functions in ( p.80|) contain Tc both 
in first and third argument so that their asymptotics at Tc — make no 
sense. 

Thus, we stop here further calculation of the Coloumb integral getting, 
however, as our main result the fact that {2'yrc)~^ should be integer number, 
in the variational approach to the model, to have finite energy of the ground 
state. We consider this as very interesting result deserving rather involved 
calculations made above. 

Also, we have a detailed technical view on the problems which arise when 
dealing with molecular integrals with the Hulten potential. Practically, this 
means that there is a very little hope that the exchange integral ( p.7| ), which 
is structurally much more complicated than the above considered Coloumb 
one, can be calculated exactly for the case of Hulten potential. 

Because of these difficulties, below we use appropriate simplified poten- 
tials, instead of Hulten potential, to have some analytical set up for the 
variational analysis of the Santilli-Shillady model. Clearly, by this we go to 
some approximation to the original Santilli-Shillady model. 

3.1.2 Coloumb integral for exponential screened Coloumb poten- 
tial 

We use simple function to mimic Hulten potential. Namely, we approximate 
the general potential ( |2.2D by 

p2 Ap-ri2/rc 

Vin,) = Vc + V, = - , (3.86) 

ri2 
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where A and rc are positive parameters. It has similar behavior both at short 
and long distances. Indeed, at long distances, ri2 — >■ oo, we can ignore Ve and 
the behavior is solely due to the Coloumb potential while its series expansion 
about the point ri2 = (short distances) is 

V^(ri2)|n.^o = + - - + 0{rl,). (3.87) 

ri2 rc 2rc 

Here, we should put A = VqTc to have the same coefficient at rfg^ in the 
ri2 — >■ asymptotics as it is in the case of Hulten potential; see Eq.(p.4|). 
Using Eq. (|3.21|) we have 

^ = ^r, = |— , /3 = 1,2,..., (3.88) 

where /3 is a number of energy levels of isoelectronium. Taking /3 = 1 we 
have, in atomic units (^ = 1, m = me/2 = 1/2), 

A = —. (3.89) 

rc 

Below, we calculate the Coloumb integral ( p.5|) , with the exponential 
screened Coloumb potential VJ= defined by Eq.( 

Tl2 ^2 



C'e = / dv,dv, |/(r„0ri/(n2)r, (3.90) 



Below, we present some details of calculation of the Coloumb integral ( |3.90| ). 
Apart from the case of Hulten potential considered in Sec. |3.1.1| , it appears 
that this integral can be calculated in terms of elementary functions. 
The integral we are calculating is 

= / dv,dv2 \firai)\Vin2)\^ (3.91) 



where 



f{r) = \r-e-^\ (3.92) 

V TT 
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and dvi and dv2 are volume elements for the first and second electron, respec- 
tively. We use spherical coordinates. In spherical coordinates {rf,2,92,^2), 
with polar axis directed along the vector ff,i, we have 



ri2 = yni^ + - '^nin2 cos 62- (3.93) 

We use these coordinates when integrating over second electron. In spherical 
coordinates (rti, 61, y^i), with polar axis directed along the vector R, we have 



Tai = V^w^ + - 2rbii?cos^i. (3.94) 

We use these coordinates when integrating over first electron. 
First, we integrate over angular coordinates of second electron, 

Ji = dip2 de2 ^e-^^^^^r^s'sin^a, (3.95) 

J J ri2 n 



where ri2 is defined by (|3.93|) . It is relatively easy to calculate this integral. 



Further, integrating on radial coordinate ri,2 must be performed in separate 
intervals, 

rbi oo 

l2 = J drh2 h{rb2 < ^6i) + j drh2 h{rb2 > rbi), (3.97) 

r-61 

where 

in, - n2y = 1'''- (3.98) 



with the result 



(3.99) 



r6i(l - 47 Vc 

Now, we turn to integrating over coordinates of first electron, (r^i, 6*1, y^i), 
J3 = y" d^^y" /2^e-2^""^rbi^sin^i, (3.100) 
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where Vai is defined by ( p.94| ). We obtain after tedious calculations 



4^ 2 

c 



R{1 - A-f\. 



2 



(3.101) 



X 



X 



Again, we must further integrate in rbi by separate intervals, 

R 00 

h = j drbi Isini <R) + j dm h{rbi > R) = hi + ^42, 

_R 

obtaining after rather tedious calculations 



(3.102) 
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(3.103) 



12i?(l - 47^^2)4 

ggg27iJ-iJ/rc^3|^(4^^27i?re+i?+re) + l)(l-27re)^+e^^-^(27r,+ l)2(47r,-l))r, 
-3e^^^ + 37(-647^(7i?(87i? + 13) + A)r^^ + 167^(7/^(247/? + 31) + 9)rc^ 

-47(7i?(247i? + 23) + Q)r^^ + i?(87/2 + 5)) - e^^-^7X 
(647^(7i?(47i?(27i? + 9) +57) +36)r,^-487=^(7i?(47i?(27i?+7) + 21) -9)rc 
+ 127(7i?(47/2(27i? + 5) + 1) - 6)rc^ - R{A-fR{2-fR + 3) - 3) + 3)) + 3 



'42 



4i?(l - 27r,)2(i + 27r,)' 



■X 



(3.104) 



X 



1-e 



47i? 



1287''i?Ve'^ + 7( (5 - 3e"^^)i? - 4(e 



4'yR\ 



l)r. 



+167^i?re^(-8i? + (16e2^^-^/"^ + 3e^^^ - 13)re) + 
+ 16e-«/^=7V/ ((86^^^ + 3e^^+«/'-= - 13e^/'^^i? - 4(e^^^ - 1) (2e2^^ - e^/'^^r,) + 
+472(2/?^ _ (3e47fi - 5)Rr, + 3(e^^^ - l)re2) + 
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In calculating I we put the condition 

67rc < 1, 



(3.105) 



which is necessary to prevent divergency at the endpoint r^i = 00. Collecting 
the above two integrals we obtain 



CL 



e-2^^(-/?(3 + 27/^(3 + 27/?)) (3.106) 



6i?(l - 47^^2)4 

+l2-i'^Rr^^{h + 27i?(5 + 27/?)) - 487^i?r/(15 + 27/^(7 + 27/?)) 



„5^ 6^-i?/r, 



+647X''(24 + 7/?(33 + 27i?(9 + 27/?)))] - 15367^6 

Thus, we have finally for the Coloumb integral for exponential screened 
Coloumb potential, 



Kl-A2)4 p 



-{p + 2p2 + p3) + 3A2(5p + 10p2 + 4p3) (3.107) 



-A^(15p + 14p2 + 4p3) +A^(8 + lip + 6p2 + -p^ 



Here, we have used notation A = 27rc, and also A < 1/3 due to Eq. (|3.105|) . 
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Figure 5: The Coloumb integral C'^ as a function of p, Eq. ( p.l07| ), at A = 
1/37. Here, p = ■jR, where R is the internuclear distance, and A = 2'jrc, 
where Vc is the correlation length parameter. 
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Figure 6: The Coloumb integral C'^ as a function of Tc, Eq. ( p.l07| ), at 
p = 1.67. For Tc > 0.2 a.u., the regularized values of C'^ are presented. 
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Figure 7: The Coloumb integral C'^ as a function of rc, Eq. ( p.l07| ), at 
p = 1.67. More detailed view. 

The total Coloumb integral is 

C'e = C'c-C (3.108) 



where Cq is wellknown Coloumb potential part given by Eq. (|3.13| ) 



Below, we turn to the other potential, Gaussian screened Coloumb po- 
tential, considered by Santilli and Shillady |jl|. The Coloumb integral for 
this potential can be calculated exactly, and the result contains one special 
function, the error function erf(z). 

3.1.3 Coloumb integral for Gaussian screened Coloumb potential 

In this Section, we calculate the Coloumb integral for the case of Gaussian 
screened potential. Namely, we approximate the general potential (|2.2|) by 



V{n2) =Vc + Vg = , (3.109) 

r 12 f 12 
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where A and c = are positive parameters. At long distances, ri2 oo, 
we can ignore Vg while its series expansion about the point ru = is 



e^-A A 



ri2 c 



+ -ri2 + 0(rf2). 



(3.110) 



Here, we should put A = VqTc to have the same coefficient at rfg^ in the 
ri2 asymptotics as it is in the case of Hulten potential; see Eq.( p.4|) . 
The Coloumb integral is 



C'q = j dvidv2 



e2 



The integral we are calculating is 



i/(r.i)ri/Mr. 



(3.111) 



C'g = j dVidV2 



ri2 



i/(rai)n/Mp 



(3.112) 



where notation and coordinate system are due to Sec. p.l.2| . First, we inte- 
grate over angular coordinates of second electron. 



27r vr 



h = J dip2 j d6- 





where ri2 is defined by (|3.93|) . We have 



2 ^e-2^^-r,2^sin52, 

ri2 vr 



(3.113) 



h = 
where 



erf(i 



erf(2;) 



erf(i 



e"* dt 



[ni - rb2) 



c 



-) , (3.114) 



(3.115) 



is error function. Further, integrating on radial coordinate rb2 must be per- 
formed in separate intervals, 

rn 00 

l2 = J drb2 h{rb2 < Hi) + J dn2 h{n2 > rbi), (3.116) 
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where 



with the result 



2^1 ni- rb2, rb2 < Tfei, 



(3.117) 



A'y^/ce 



-2-yrti-rti^/c 



4r, 



47v^( 



bl 



(3.118) 



:i + 27(rw - C7))(erfc(^^^^^) + 2erfc(v^7) - 2) 



+e^^^" (27(rbi + 07) - l)erfc( 



ni + C7 . 



where erfc(2;) = 1 — erf(z). Now, we turn to integrating over coordinates of 
first electron, {rbi,Oi,{pi), 



n 2n 



h = J dOiJ dipi /2— e-'^'-'-Vw^sin^i, (3.119) 


where r^i is defined by ( p.94| ). We obtain after tedious calculations 

T _ ^Vc'T^„-IiL-27(v/(fl-r6i)2+2r6i + V(«+'-M)2) 



(3.120) 



X 



+27e2^(V(^-'-")'+'^"y (i? + my 



^e^+^^' ( (1 + 27rbi - 2c72)(2erfc(7v^) + erfc(-'''' ^""^ 



-2) 



+ (1 + 27r6i + 2c7^)e^^'^"erfc( 



Again, we must further integrate in r^i by separate intervals, 

R oo 



h = j dm him <R)+J dm hini > R). (3.121) 
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First, we replace the endpoint Vhi = oo by finite value Vhi = A to avoid 
divergencies at intermediate calculations. After straightforward but tedious 
calculations we obtain rather long expression so that we do not represent it 
here noting however that the following integrals are used during the calcula- 
tions: 

/ eri(z)dz = ^ + z erf(z), (3.122) 
J ^/tt 

[ z eTf{z)dz = - ^erf(z) + J-z^erf(z), (3.123) 

J 2y7r 4 2 

/ e-'''eri{z)dz = --e-"^erf(^) + -e"'/^erf(^ + z), (3.124) 

J Qj Oj Zi 

I ze-'''ei:i(z)dz = ^e""^-^' - ^e-^Ul + az)eri(z) (3.125) 

1 V-l)e"'/Vf(^ + ^), 



2a2^ ' '2 



Je-^^-'^'dz = ^e'^'/(^'')erf(^±|^), (3.126) 

ze-^'-'-'dz = -^e— - ^e«'/(^^)erf(^^±|^). (3.127) 

Using limA_>oo erf(A) = 1 and replacing welldefined exponentially decreasing 
terms by zero, we obtain some finite terms and big number (about fourty) of 
A dependent terms, which are unbounded at A — > oo. All the divergent terms 
totally cancel each other so the final expression turns out to be automatically 
finite. 

As the result, we obtain the Coloumb integral for Gaussian screened 
Coloumb potential in the following form: 



^ 96p 



-(60 + 96p + 48p^)K + (32 + 48p)«;'^ - 16«;^ (3.128) 



+ ((60 + 16p^)«; - 32k^ + 16K^)e^''"^ 
+V^e''' ( (30p + 8p3 - 36pK^ + 24p«;^) (2erf («;) - erfc(^ - k) - e^^erfc(^ + k)) 

\ K K 
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Figure 8: The Coloumb integral Cg as a function of p, Eq. ( |3.128| ), at 2k = 
A = 1/37. Here, p = 'jR, where R is the internuclear distance, and A = 2'yrc, 
where is the correlation length parameter. 



+ (15+24p2-(18+24p2)/t2+12/t^-8K'^)(2erf(/t)-erfc(--/t)+e^''erfc(-+/€)) 



K 



where we have used notation 



K = 7VC = 



(3.129) 



The total Coloumb integral is 



C — C ~ C 



(3.130) 



where C'q is given by Eq. (p.l3|) 
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Figure 9: The Coloumb integral Cg as a function of r^, Eq. ( p.l28| ), at 
p = 1.67. 
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Figure 10: The Coloumb integral C'g as a function of r^, Eq. ( |3.128| ), at 
p = 1.67. More detailed view. 
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3.1.4 Exchange integral 

Our general remark is that all calculations for the above Coloumb integrals 
are made in spherical coordinates, which correspond to spherical symme- 
try of the charge distributions of both Is electrons, iV'l^ai)!^ and \4'{ra2)\'^, 
each moving around one nucleus. One can use prolate spheroidal coordi- 
nates, which are exploited sometimes when integrating over coordinates of 
last electron, but we have encountered the same problem of big number of 
terms in the intermediate expressions, with no advantage in comparison to 
the use of spherical coordinates. 

Unlike to Coloumb integral, calculation of exchange integral should be 
made in the spheroidal coordinates, which correspond to spheroidal symme- 
try of charge distributions of the electrons, 4'*{rai)ip{ri,i) and ip*{ra2)4'{rh2), 
each moving around both the nuclei, a and b. 

Calculation of the exchange integral. 



essentially depends on the form of the potential V{ri2) in the sense that 
the integration can be made only in s;)/ieroz(ia/ coordinates, {xi,yi,(pi) and 
(a;2, 2/2, ¥^2), and one should use an expansion of V{ri2) in the associated 
Legendre polynomials. 

For the usual Coloumb potential, V{ri2) = rf2^, it is rather long (about 
12 pages to present the main details) and nontrivial calculation, where Neu- 
mann expansion in terms of associated Legendre polynomials, in spheroidal 
coordinates, is used (celebrated result by Sugiura, see Eq.( |3.15| )). 

In general, any analytical square integrable function can be expanded in 
associated Legendre polynomials. However, in direct calculating of the ex- 
pansion coefficients by means of integral of the function with Legendre poly- 
nomials, one meets serious problems even for simple functions. Practically, 
one uses, instead, properties of special functions to derive such expansions. 

We mention that there is Gegenbauer expansion 0, having in a particular 
case the form P 



where ji{z) and ni (z) are spherical Bessel and spherical Hankel functions of 
first kind, respectively, 612 is angle between vectors fi and r2, and ri = |ri|. 




(3.131) 




(3.132) 
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f2 = 1^2!; ^1 < f2- Spherical harmonics Yifi{6i2) can be rewritten in terms of 
Legendre polynomials due to the summation theorem. 

We note that this expansion can be used, at A; = i/rc, to reproduce ex- 
ponential screened potential, Ve(ri2), and to calculate associated exchange 
integral (p.l31|) but, however, it concerns spherical (not spheroidal) coordi- 
nates, (ri,^i,v2i) and (r2, 6^2, ^^'2)- 

For Hulten potential V/j(ri2), exponential screened potential, Ve{ri2), and 
Gaussian screened potential, V^(ri2), which are of interest in this paper, we 
have no such an expansion in spheroidal coordinates. To stress that this 
is not only the problem of changing coordinate system, we mention that 
the solution of usual 3-dimensional wave equation, Aip + k'^ip = 0, is given 
by function e*'^^/r, in spherical coordinates, to which one can apply Gegen- 
bauer expansion, while in spheroidal coordinates its solution is represented 
by complicated function containing infinite series of recurrent coefficients [0]; 
see also Sec. 3.4. As the result, we have no possibility to calculate exactly 
exchange integrals for these non-Coloumb potentials. 

In order to obtain approximate expression for the exchange integral for the 
case of the above non-Coloumb potentials, we make analysis of asymptotics 
of the standard exchange integral (i.e. that for the Coloumb potential), 
Eq. (|3.15|) . It is easy to derive that 

^qp^oo-e-^ (3.133) 
at long distances between the nuclei, and 

^qp=o = ^7, (3.134) 

in the case of coinciding nuclei. At 0, we should have the same 

asymptotics for exchange integral for each of the above non-Coloumb poten- 
tials because these potentials behave as Coloumb potential at 0. 

In both the limiting cases, p ^ 00 and p = 0, the exchange integral for the 
non-Coloumb potentials is simplified, and one can use spherical coordinates 
since the two-center problem is reduced to one-center problem. We consider 
two limiting cases. 

a) p = 00 case. 

This case is trivial because exchange integral tends to zero due to lack of 
overlapping of the wave functions of two H atoms. 
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b) p = case. 

In this case , we have Vai = r^i = ri and = rb2 = T2 so that Eq.( p.l3l| ) 
becomes 

S' = I dv^dv2 V(ri2)|/(ri)n/(^2)|', (3.135) 

One can see that this is the case of He atom with two electrons in the ground 
state. Evidently, in terms of our anzatz ( |3.1| ) we have complete overlapping 
of the wave functions. 

Even the above mentioned simplification of the exchange integral and use 
of spherical coordinates does not enable us to calculate straightforwardly the 
integral (|3.135|) for the non-Coloumb potentials, Vh, Ve, or Vg] the integrands 
are still too complicated. This indicates that we should use expansion of these 
potentials in Legendre polynomials, in spherical coordinates, to perform the 
integrals. Only exponential screened potential is given such an expansion 
here. Namely, this is Gegenbauer expansion (|3.132|) , owing to which we can 
calculate the exchange potential for exponential screened potential V^, to 
which we turn below. 



Exchange integral for the exponential screened Coloumb poten- 
tial Ve, at p = 0. 



The integral is 



|p=0 



5 f 

-7 - y dvidv2 



ri2 



l/(ri)ri/(r2)P, (3.136) 



where we have used Eq. (|3.134|) for the usual Coloumb potential part of the 
integral. In the Gegenbauer expansion (p.l32 ), we assume k = i/rc to repro- 
duce the potential Ve^ru). Since the wave functions /(ri) and /(r2) given by 
Eq.( |3.10|) do not depend on the angles, only / = 0, m = term of the expan- 
sion (|3.132|) contributes to the exchange integral (|3.136|) due to orthogonality 
of Legendre polynomials. Using 



sm z, 



-le 




we thus have 



ri2 



Y^— sin krie^^^'^ 
sin fcroe*'^''! 

krrr2 



0,0 



ri < r2, 
ri > r2. 



(3.137) 



(3.138) 
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Then the exchange integral ( |3.136| ) is written as 



F' 



S|p=0 



A 



-7- U7ir^dr2 UTrrfdri- smkne"'''^—e-'^''—e 



Jkr2 -2'yri -2'yr2 



krir2 



IT 



IT 



(3.139) 



A 



+ UTrrfdri- sin A;r2e*'="i^e-^^"i^e 



ikri —2'yri 7 ^—2'yr2 



krir2 



TT 



where Anrf and 47rr2 are volume factors. The two above integrals over ri can 
be easily calculated, with the result 



+ 47 



and 



2^2 



k 

— (fc^ + 47^) (/c COS kr2 + 27 sin kr2 



■^1 + (27 — iA;) sin/cr2e* 



Summing up these terms and integrating over r2 we get after some algebra 

^ A73 



F' 



-e;|p=o 



:7 + 



2{k + 2^7)' 



;(r + 8iA;7-207^). (3.142) 



Inserting 



(3.143) 



to reproduce the potential and denoting A = 27rc we write down our final 
result, 



F' 



o7 



7AA2 



Note that, at 0, i.e. at A ^ 00, we have 



1 + 4A + 5A2). 



(3.144) 



F' 



E\p={i 



-7 A'-j 



(3.145) 



that is in agreement with the value ( p.l34| ). We should to emphasize here 
that Eq.( p.l4'^ ) is exact result for the exchange integral £'^, at p = 0. 
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Figure 11: The exchange integral as a function of p, Eq. ( p.l46| ), at 
A = 1/37. Here, p = 'jR, where R is the internuclear distance, and A = 2'yrc, 
where is the correlation length parameter. 



Next step is to implement p dependence into ( |3.144| ) following to natural 
criteria. To restore partially p dependence in the exchange integral (|3.144| ), 
we use exact result ( |3.15| ), and write down for the p dependent exchange 
integral the following approximate expression: 

g' g' g' r^g' — + -A + -A^)-^' , (3.146) 

where S'^ is standard exact exchange integral for Coloumb potential given by 
Eq. (|3.15|) while the approximate A dependent part arised from our potential 
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Figure 12: The exchange integral 8^ as a function of Tc, Eq. ( 3.146 ), at 
p = 1.67. 
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Figure 13: The exchange integral 8'^ as a function of r^, Eq. ( |3.146| ), at 
p = 1.67. More detailed view. 
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We have a good accuracy of the approximation ( p.l46| ). Indeed, exchange 
integrals make sensible contribution to the total molecular energy at deep 
overlapping of the wave functions, S > 0.5, and we have made calculation 
just for the case of complete overlapping, S = 1, with the necessary asymp- 
totic factor, e"^'', provided by E'cip)- Note that at A — oo, the term 8'^ of 
Eq. (|3.146| ) becomes AE'^^, as it should be because at A ^ oo (no screening) 



we have Ajryi- In addition, although there is no possibility to restore 

completely p dependence for the second term in r.h.s. of Eq. (|3.146|) , we have 
got information on A dependence, which is of most interest here. 

3.2 Numerical calculations for the Vg-based model 

In this Section, we consider the case of exponential screened potential Ve = 
^g-''i2/rc^^^ for which we have calculated all the needed molecular integrals. 
The H2 molecule energy, due to Eq.(|3.2|), is written as 



F i-^: R Ar)-2 ^^^'^ 2{C + SS) - {C'g ~ + S'g - S'^) 1 
E_,(7, R, A, r,) - 2^^^^ + -, 

(3.147) 

where the specific terms are the Coloumb integral C'^ given by Eq.( ^.107| ) and 
the exchange integral given by ( p.l46|) . We should find extremum of Emoi 



as a function of our basic parameters, 7, R, A, and r^. We are using notation 
p = '-fR and A = 2'~frc so that our four parameters are 7, p, A, and A. In 
general, the number of energy levels of isoelectronium can also be viewed 
as a parameter of the model. However, we restrict our consideration by the 
one-level case, = 1; see Sec. |2?1. 



3.2.1 Minimization of the energy 



First, we analyze the A dependence of Emoi- Due to Eq. (|3.89|) , for one-level 



isoelectronium we have A = r^^, that can be identically rewritten as 

A = ^. (3.148) 
A 

Thus the A dependence converts to 7 and A dependence. This is the conse- 
quence of consideration of the Hulten potential interaction for the electron 



pair made in Sec. 2.1 
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Second, we turn to 7 dependence. Due to ( p.l48| ), the A dependent parts, 
Cg and S'^, acquire additional 7 factor and thus become 7^ dependent. The 
other molecular integrals depend on 7 linearly so that we define accordingly, 

r - -C f - -f C" - -C F' - - f C" - —C F' - —F' 

(3.149) 

Inserting the computed integrals A and A' into ( p.l47| ) we have 

Emoih,p,X) = -a^ + b^^, (3.150) 

where 

2 + 2C + 4SS-C'c-S'c 1 
a(p,A) = ^^-^^ (3.151) 

and 

Kp. A) = ^^^2 ■ (3.152) 

The value of 7 corresponding to an extremum of Emoi is found from the 
equation dEmoi/d'j = 0, which gives the optimal value 

7opt = ^- (3.153) 



Inserting this into ( |3.150| ) we get the extremal value of E., 



mol ) 



E^oi{p,X) = -^. (3.154) 

Using definitions of a and b we have explicitly 

_ l-2p + S^ + p{-2C-ASS + C'c + C'a) 

2p{-l + S^-2SS^ + C', + S'J ^ ^ 

and 

Next, we turn to the extremum in the parameter p. The p dependence, 
as well as the A dependence, of Emoi is essentially nonalgebraic so that we 
are forced to use numerical calculations. 
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Figure 14: The total energy E = E^oi as a function of p, Eq. ( p.l56| ), at 
A = 1/60,1/40,1/20,1/10,1/5. The lowest plot corresponds to A = 1/5 
(p = 7i?, A = 27r,). 

It appears that the A dependence does not reveal any local energy mini- 
mum while the p dependence does. Below, we use the condition, A~^ = in- 
teger number, obtained during the calculation of the Coloumb integral with 
Hulten potential Vh] see Eq. ( p.57| ). Although there is obviously no necessity 
to keep this condition for the case of exponential screened potential Ve, we 
consider it as a prescription for allowed values of A. 

Since the A dependence of the energy has no minimum we can use fitting 
of the predicted energy EmoiW to the experimental value by varying A. This 
allows us to estimate the value of the parameter A, and thus the value of the 
effective radius of the isoelectronium = X/2'-fopt- 

3.2.2 Fitting of the energy and the bond length 

The procedure is the following. We fix some numerical value of A, and iden- 
tify minimal value of Emoi{p,X), given by Eg. ( p. 15^ ), in respect with the 
parameter p. This gives us minimal energy and corresponding optimal value 
of p, at some fixed value of A. Then, we calculate jopt by using Eq. (|3.155|) , 
and use obtained values of popt and 'jopt to calculate values of Ropt and Tc. 
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We calculated minimal values of Emoi in P, for a wide range of integer 
values of A~^. The results are presented in Tables 2 and 3, and Figures 15 
and 16. One can see that the energy Emoi decreases with the increase of Tc 
(proportional to size of isoelectronium), as it was expected to be. 

We note that all the presented values of Emoi in Tables 2 and 3 are lower 
than that, E^^j^ = —1.139 a.u., obtained via two-parametric Ritz variational 
approach to the standard model of H2 (see, e.g., [H), which is the model 
without the assumption of short-range attractive potential between the elec- 
trons. This means that the VJ=-based model gives better prediction than the 
one of the standard model, for any admitted value of the effective radius of 
isoelectronium Vc > 0. Indeed, the standard prediction -E^JJ = —1.139 a.u. 
is much higher than the experimental value Eexper[H2] = —1.174474 a.u. 

3.2.3 The results of fitting 
Best fit of the energy Emoi- 

Due to Table 2 (see also Fig. 15), the experimental value, Eexp[H2] = 
— 1.174... — 1.164 a.u. (here we take 0.9% uncertainty of the experimental 
value) is fitted by 

= 0.0833.. .0.0600 a.u., (3.157) 

i.e. A = 1/5. ..1/7, with the optimal distance, Ropt = 1.3184. ..1.3441 a.u. 
We see that the predicted Ropt appeared to be about 6% less than the ex- 
perimental value Rexper[H2] = 1.4011 a.u. We assign this discrepancy to the 
approximation we have made for the exchange integral (|3.146|) . 
Below we fit Ropt, to estimate the associated minimal energy. 

Best fit of the internuclear distance R. 

Due to Table 2 (see also Fig. 16), the experimental value of the inter- 
nuclear distance, Rexp = 1.4011 a.u., is fitted by Vc = 0.0115 a.u., with the 
corresponding minimal energy E^in = —1.144 a.u., which is about 3% big- 
ger than the experimental value. Again, we assign this discrepancy to the 
approximation we have made for the exchange integral ( |3.146D , and take 

= 0.0115 a.u., (3.158) 
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i.e. A = 1/37, as the result of our final fit noting that (a) in ref. the value 
Tc = 0.0112 a.u. has been used to make exact numerical fit of the energy, 
with corresponding R= 1.40 a.u., and (b) we have less discrepancy. 

The weight of the pure isoelectronium phase. 

To estimate the weight of the pure isoelectronium phase, which can be 
viewed as a measure of stability of the pure isoelectronium state, we use the 
above obtained fits and the fact that this phase makes contribution to the 
total molecular energy via the Coloumb and exchange integrals. 

According to Eg. ( p.l47| ), the isoelectronium phase displays itself only by 
the term Pe = |Ce(7,p, A) +^g(7,p, A)| while the Coloumb phase displays 
itself by the corresponding term Pc = \C'(j{j,p) + £^^(7,p)|. Putting the 
total sum Pc + Pe = 1, i-e. Pc + Pe is 100%, the weights are defined simply 
by 

The weight for the best fit of R. 
At the values A = 1/37 (i.e. r, = 0.0115), 7 = 1.1706, and p = 1.6320, for 
which we have minimal Emoi = —1.144 and optimal R = 1.40, we get the 
numerical values of the weights. 

We = 0.84% (3.160) 

for the pure isoelectronium phase, and Wc = 100% — = 99.16% for the 
Coloumb phase. 

The weight for the best fit of Emoi ■ 
At the values A = 1/5 (i.e. = 0.0833 a.u.), 7 = 1.2005, and p = 1.5827, 
for which we have minimal Emoi = —1.173 a.u. and optimal R = 1.318 a.u., 
we obtain 

We = 6.16%, Wc = 93.84%. (3.161) 

From the above two cases, one can see that the weight of pure isoelectronium 
phase is estimated to be 

We ~ 1...6%, (3.162) 

for the predicted variational energy Emoi = —1.143... — 1.173 a.u. 

The biggest possible weight. 
Note that in our V^-based model the biggest allowed value of A is A = 1/4 
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(i.e. Tc = 0.1034) because A < 1/3, to avoid divergency of the Coloumb 
integral Ce- For this value of A, we obtain minimal Emoi = —1.182 a.u. and 
optimal R = 1.297 a.u. This value corresponds to the ^i^'^'est possible weight 
of the pure isoelectronium phase, 

We = 7.32%, (3.163) 

within our approximate model. 

The following three remarks are in order. 

(i) We consider the existence of this upper limit. We < 7.32%, as a highly 
remarkable implication of our V^-model noting however that it may be artifact 
of the use of the exponential screened Coloumb potential. 

(ii) Another remarkable implication is due to the condition, A~^ = integer 
number, obtained for the case of Hulten potential. One can see from Table 2 
that the energy Emoi varies discretely with the discrete variation of A~^. This 
means that there is no possibility to make a "smooth fit". For example, at 
A = 1/5, we have Emoi = —1-173, and the nearest two values, A = 1/4 and 
A = 1/6, give us Emoi = —1.182 and Emoi = —1.167, respectively. Therefore, 
owing to the above condition the model becomes more predicitive. 

(iii) Numerical calculation shows that the formal use of the exact Coloumb 
integral C'g, given by Eq.( |3.12g| ), of the Gaussmn screened Coloumb potential, 



instead of Cg, in Eq.( p.l47D gives us approximately the same fits. Namely, 
the best fit of the energy is achieved at A = 1/5, with rc = 0.1042, optimal 
R = 1.323, and minimal E^oi = -1.172. Also, the best fit of = 1.40 
is at A = 1/29, for which rc = 0.0147 and minimal Emoi = —1.144. Here, 
we have used the same exchange integral as it is for the case of exponential 
screened potential so these fits have been presented just for a comparison with 
our basic fits, and to check the results. Note that for the case of Gaussian 
screened Coloumb integral we have no restriction on the allowed values of A. 
Analysis shows that, at big values of A, e.g. at A > 4, the integral Cg, given 
by (|3.128|) , rapidly oscillates in the region of small p {p < 0.5). This means 
that when the correlation length rc becomes comparable to the internuclear 
distance an effect of instability of the molecule arises. This can be viewed as 
a natural criterium to fix the upper limit of A. Normally, we use the values 
A < 1, for which case there are no any oscillations of Cg (see Fig. 9). 
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\ —1 


Tc, a.u. 




■^mini a-U. 


4 





10337035071618050 


1 


297162129235449 


-1 


181516949656805 


5 





08329699109108888 


1 


318393698326879 


-1 


172984902150024 


6 





06975270534273319 


1 


333205576478603 


-1 


167271240301846 


7 





05999677404817234 


1 


344092354783681 


-1 


163188554065554 


8 





05263465942162049 


1 


352417789644028 


-1 


160130284706318 


9 





04688158804756491 


1 


358984317233049 


-1 


157755960428922 


10 





04226204990365446 


1 


364292909163710 


-1 


155860292450436 


11 





03847110142927672 


1 


368671725082009 


-1 


154312372623724 


12 





03530417706681329 


1 


372344384866235 


-1 


153024886026671 


13 





03261892720535206 


1 


375468373051375 


-1 


151937408039373 


14 





03031323689615631 


1 


378157728092548 


-1 


151006817317425 


15 





02831194904031777 


1 


380497017045902 


-1 


150201529091051 


16 





02655851947236431 


1 


382550255552670 


-1 


149497886394651 


17 





02500959113834722 


1 


384366780045693 


-1 


148877823925501 


18 





02363136168905809 


1 


385985219224291 


-1 


148327310762828 


19 





02239708901865092 


1 


387436244558651 


-1 


147835285349041 


20 





02128533948435381 


1 


388744515712491 


-1 


147392910500336 


21 





02027873303335994 


1 


389930082626193 


-1 


146993041730378 


22 





01936302821907175 


1 


391009413196452 


-1 


146629840949675 


23 





01852644434336641 


1 


391996158084790 


-1 


146298491232105 


24 





01775915199935013 


1 


392901727808297 


-1 


145994983116511 


25 





01705288514774330 


1 


393735733699196 


-1 


145715952370148 


26 





01640064219648127 


1 


394506328745493 


-1 


145458555325045 


27 





01579645313764336 


1 


395220473843219 


-1 


145220372020229 


28 





01523519631632570 


1 


395884147817973 


-1 


144999330178493 


29 





01471245291356761 


1 


396502514589167 


-1 


144793644973560 


30 





01422439038752817 


1 


397080057337240 


-1 


144601770891686 



Table 2: The total minimal energy Emin and the optimal internuclear distance 
Ropt as functions of the correlation length r^. The exponential screened 
Coloumb potential case (see Figures 15 and 16). 
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A ^ 


rc, a.u. 




■^mini S-U. 


31 





01376766836566138 


1 


397620687025853 


-1 


144422362947838 


32 





01333936209977966 


1 


398127830817745 


-1 


144254245203342 


33 





01293689977547854 


1 


398604504597664 


-1 


144096385030938 


34 





01255801083612469 


1 


399053372836414 


-1 


143947871939897 


35 





01220068312791624 


1 


399476798299823 


-1 


143807900045981 


36 





01186312715793131 


1 


399876883556063 


-1 


143675753475045 


37 





01154374612489787 


1 


400255505817128 


-1 


143550794143290 


39 





01095393745919852 


1 


400954915288619 


-1 


143320213707519 


40 





01068107105944273 


1 


401278573036792 


-1 


143213620508321 


41 





01042146833640030 


1 


401586548200467 


-1 


143112256673494 


42 





01017418516195214 


1 


401879953246168 


-1 


143015746732479 


43 





00993836493541500 


1 


402159797887369 


-1 


142923750307661 


44 





00971322867044429 


1 


402427000676349 


-1 


142835958109381 


45 





00949806639934841 


1 


402682399061957 


-1 


142752088467028 


46 





00929222969498477 


1 


402926758144872 


-1 


142671884314343 


47 





00909512514431396 


1 


403160778323019 


-1 


142595110561057 







00890620863525624 


1 


403385101987775 


-1 


142521551794315 


49 





00872498034101540 


1 


403600319405678 


-1 


142451010262626 


50 





00855098030451296 


1 


403806973898863 


-1 


142383304102633 


51 





00838378454080327 


1 


404005566419838 


-1 


142318265775268 


52 





00822300158793934 


1 


404196559601683 


-1 


142255740683024 


53 





00806826944722482 


1 


404380381352424 


-1 


142195585944305 


54 





00791925286251402 


1 


404557428052374 


-1 


142137669304475 


55 





00777564089552400 


1 


404728067404676 


-1 


142081868166104 


56 





00763714476025456 


1 


404892640982100 


-1 


142028068723488 


57 





00750349588477794 


1 


405051466507240 


-1 


141976165188595 


58 





00737444417302681 


1 


405204839898059 


-1 


141926059097351 


59 





00724975644291090 


1 


405353037106507 


-1 


141877658686723 


60 





00712921502024112 


1 


405496315774223 


-1 


141830878334298 



Table 3: The total minimal energy Emin and the optimal internuclear distance 
Ro-pt as functions of the correlation length r^- The exponential screened 
Coloumb potential 14 case (see Figures 15 and 16). 
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Figure 15: The total minimal energy E^^n as a function of the correlation 
length Tc- The exponential screened Coloumb potential Ve case (see Tables 2 
and 3). 




Figure 16: The optimal internuclear distance Rgpt as a function of the cor- 
relation length Tc- The exponential screened Coloumb potential Ve case (see 
Tables 2 and 3). 
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